Load data

raw <- 
  bind_rows(
    read_excel(file.path("data", "2017_0413_0822_GB_data_commpilation.xlsx")) 
  ) %>% 
  select(File, Nr., Start, Rt, End, # filename and retention times
         Ampl44=`Ampl 44`, Ampl45=`Ampl 45`, # amplitudes
         Area44=`rIntensity 44`, Area45=`rIntensity 45`, # areas
         Intensity=`Intensity All`, #peak intensities
         R13C = `R 13C/12C`, d13C = `13C/12C` # ratio and delta values
  )

Map peaks

### CHANGE MAPPING FILE NAME ###
data_matched <- map_peaks(raw, metadata_file = file.path("metadata", "smoky_hollow_peaks_April13_August22_compiled.xlsx"), quiet = FALSE)
## INFO: 1010 peaks in 50 files successfully assigned.
## WARNING: 13 file(s) have no (or undefined) peak map assignments
## 1668__170820_A5_PTV_3_OG038 Ad.dxf
## 1669__170820_A5_PTV_3_OG052 Ad.dxf
## 1671__170820_A5_PTV_3_OG044 Ad.dxf
## 1672__170820_A5_PTV_3_OG014 Ad.dxf
## 1674__170820_A5_PTV_3_OG041 Ad.dxf
## 1675__170820_A5_PTV_3_OG058 Ad.dxf
## 1678__170820_A5_PTV_3_OG061 Ad.dxf
## 1679__170820_A5_PTV_3_OG042 Ad.dxf
## 1681__170820_A5_PTV_3_OG046 Ad.dxf
## 1682__170820_A5_PTV_3_OG051 Ad.dxf
## 1684__170821_A5_PTV_3_OG039 Ad.dxf
## 1685__170821_A5_PTV_3_OG043 Ad.dxf
## 1686__170821_A5_PTV_3_OG042 Ad.dxf
## WARNING: 610 peaks are unassigned or missing
## # A tibble: 610 × 22
##                                                               File   Nr.
##                                                              <chr> <dbl>
## 1  1666__170820_A5_PTV_0_5_A5dilutionsTest_CO2_newVial_14ng@ul.dxf    NA
## 2  1666__170820_A5_PTV_0_5_A5dilutionsTest_CO2_newVial_14ng@ul.dxf    NA
## 3  1666__170820_A5_PTV_0_5_A5dilutionsTest_CO2_newVial_14ng@ul.dxf    NA
## 4  1666__170820_A5_PTV_0_5_A5dilutionsTest_CO2_newVial_14ng@ul.dxf    NA
## 5  1666__170820_A5_PTV_0_5_A5dilutionsTest_CO2_newVial_14ng@ul.dxf    NA
## 6  1666__170820_A5_PTV_0_5_A5dilutionsTest_CO2_newVial_14ng@ul.dxf    NA
## 7  1666__170820_A5_PTV_0_5_A5dilutionsTest_CO2_newVial_14ng@ul.dxf    NA
## 8  1666__170820_A5_PTV_0_5_A5dilutionsTest_CO2_newVial_14ng@ul.dxf    NA
## 9  1666__170820_A5_PTV_0_5_A5dilutionsTest_CO2_newVial_14ng@ul.dxf    NA
## 10 1666__170820_A5_PTV_0_5_A5dilutionsTest_CO2_newVial_14ng@ul.dxf    NA
## # ... with 600 more rows, and 20 more variables: Start <dbl>, Rt <dbl>,
## #   End <dbl>, Ampl44 <dbl>, Ampl45 <dbl>, Area44 <dbl>, Area45 <dbl>,
## #   Intensity <dbl>, R13C <dbl>, d13C <dbl>, compound <chr>,
## #   is_ref_peak <chr>, map <chr>, Rt_target <dbl>, type <chr>,
## #   process <chr>, notes <chr>, depth <dbl>, n_matches <dbl>,
## #   n_overlapping <dbl>

Add standard values

standards <- read_excel(file.path("metadata", "gc_irms_indiana_A6.xlsx"))
data_w_stds <- data_matched %>% 
  filter(type == "standard", is_ref_peak == "no") %>% 
  left_join(standards, by = "compound") %>% 
  mutate(is_std = !is.na(true.d13C) | !is.na(true.d2H))

Evaluate standards

Visualize

data_w_stds %>% 
  ggplot() +
  aes(x = true.d13C, y = d13C, color = File) + 
  geom_smooth(method = "lm", se = FALSE, alpha = 0.5) +
  geom_point() +
  theme_bw() +
  theme(legend.position = "none")

Overview

std_corrections <- 
  data_w_stds %>% 
  group_by(File) %>% 
  do({
    m <- lm(d13C ~ true.d13C, data = .)
    data_frame(
      offset = coefficients(m)[1],
      slope = coefficients(m)[2],
      delta_ref_CO = -offset/slope,
      max_residual = max(summary(m)$residuals),
      rss = sum(summary(m)$residuals^2),
      r2 = summary(m)$r.squared
    )
  }) 
std_corrections %>% knitr::kable(d = 3)
File offset slope delta_ref_CO max_residual rss r2
1670__170820_A5_PTV_1_0_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 11.415 1.019 -11.204 0.300 0.498 0.994
1673__170820_A5_PTV_1_5_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 11.731 1.025 -11.440 0.241 0.279 0.997
1676__170820_A5_PTV_1_0_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 11.755 1.025 -11.463 0.347 0.408 0.995
1677__170820_A5_PTV_1_5_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 11.293 1.012 -11.155 0.230 0.252 0.997
1680__170820_A5_PTV_1_5_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 12.557 1.055 -11.904 0.206 1.058 0.988
1683__170820_A5_PTV_1_0_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 11.058 1.005 -11.004 0.345 0.425 0.994
1687__170821_A5_PTV_0_5_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 12.884 1.066 -12.083 0.377 2.747 0.969
1688__170821_A5_PTV_1_5_A5dilutionsTest_CO2_newVial_14ng@ul.dxf 11.130 1.006 -11.064 0.247 0.374 0.995
473__A5__1_5UL_30ng.dxf 10.279 1.011 -10.170 0.422 0.428 0.994
474__A5__1_5UL_30ng.dxf 10.305 1.013 -10.168 0.297 0.483 0.994
479__A5__1_5UL_30ng.dxf 10.436 1.015 -10.279 0.282 0.303 0.996
484__A5__1_5UL_30ng.dxf 9.976 0.997 -10.006 0.324 0.286 0.996

isotope values check by rentention time

data_w_stds %>% 
  ggplot() +
  aes(x = Rt, y = d13C, color = File) + 
  geom_point() +
  theme_bw() +
  theme(legend.position = "none") 

ggplotly(ggplot2::last_plot() + theme(legend.position = "none"))

peak intensities check by rentention time

data_w_stds %>% 
  ggplot() +
  aes(x = Rt, y = Intensity, color = File) + 
  geom_point() +
  theme_bw() +
  theme(legend.position = "none") 

ggplotly(ggplot2::last_plot() + theme(legend.position = "none"))

Data

message(sprintf("Infered reference tank isotopic composition: %.2f +/- %.2f",
                std_corrections$delta_ref_CO %>% mean(),
                std_corrections$delta_ref_CO %>% sd()))
## Infered reference tank isotopic composition: -10.99 +/- 0.70
# simplifed offset=ref gas correction (TODO: consider signal strengh (linearity correct), consider drift and evaluate standards between samples rather than average)
offset_avg <- std_corrections$offset %>% mean()
slope_avg <- std_corrections$slope %>% mean()

All

data_matched %>% 
  filter(type == "sample", is_ref_peak == "no") %>% 
  mutate(d13C.corrected = (d13C - offset_avg)/(slope_avg)) %>% 
  ggplot() +
  aes(compound, y = d13C.corrected, color = File) +
  geom_point() +
  theme(axis.text.x = element_text(angle = 90))

ggplotly(ggplot2::last_plot() + theme(legend.position = "none"))

Against depth

data_matched %>% 
  filter(type == "sample", is_ref_peak == "no") %>% 
  mutate(d13C.corrected = (d13C - offset_avg)/(slope_avg)) %>% 
  ggplot() +
  aes(x = depth, y = d13C.corrected, color = File, size = Area44) +
  geom_point() +
  scale_x_reverse() + 
  facet_wrap(~compound, scales = "free") + 
  coord_flip()

ggplotly(ggplot2::last_plot() + theme(legend.position = "none"))